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Abstract 

Hydrodynamic modelling is an important tool for the development of tidal 
stream energy projects. Many hydrodynamic models incorporate the effect of 
tidal turbines through an enhanced bottom drag. In this paper we show that 
although for coarse grid resolutions (kilometre scale) the resulting force exerted 
on the flow agrees well with the theoretical value, the force starts decreasing 
with decreasing grid sizes when these become smaller than the length scale of 
the wake recovery. This is because the assumption that the upstream velocity 
can be approximated by the local model velocity, is no longer valid. Using linear 
momentum actuator disc theory however, we derive a relationship between these 
two velocities and formulate a correction to the enhanced bottom drag formula¬ 
tion that consistently applies a force that remains closed to the theoretical value, 
for all grid sizes down to the turbine scale. In addition, a better understanding 
of the relation between the model, upstream, and actual turbine velocity, as 
predicted by actuator disc theory, leads to an improved estimate of the usefully 
extractable energy. We show how the corrections can be applied (demonstrated 
here for the models MIKE 21 and Fluidity) by a simple modification of the drag 
coefficient. 


* Corresponding author 

Email address: s.kramer@imperial.ac.uk (Stephan C. Kramer) 


Preprint submitted to Renewable Energy 


June 12, 2015 



Keywords: tidal turbines, tidal power, tidal stream, enhanced bottom drag, 
hydrodynamic modelling, energy resource assessment 


1. Introduction 

One of the key advantages of tidal energy as a renewable energy source, 
is the predictable nature of the resource. Methods for the detailed prediction 
of tidal dynamics using hydrodynamic numerical models have been developed 
over many years and have been applied for many different purposes. Less well 
understood is how the placement of tidal energy converters in the flow will 
modify the existing tidal currents at both local and regional scales [T]. The 
challenge here is that the detailed flow around a turbine is a three-dimensional 
phenomenon comprising far smaller length scales than those of the underlying 
tidal resource. A typical approach therefore is to model the turbine scale flow 
in a three-dimensional CFD simulation based on a actuator disc, blade element, 
or actuator-line model (see e.g. Sun et al. (22], Harrison et al. [TO], Batten 
et al. [2], Malki et al. (T5|, Churchfield et al. [3]). The effects of the turbine in 
a large scale hydrodynamic model are then parameterised, based on properties 
extracted from the CFD model. 

The main property of the turbine that needs to be parameterised is the 
amount of thrust force exerted by the turbine on the flow (and vice-versa) as 
a function of the flow speed. This also determines the amount of energy taken 
out of the flow. Thrust curves typically take the form of a quadratic function 
of current speed with a non-dimensional thrust coefficient, and can be derived 
as described above in a high-resolution CFD model, or in lab experiments. 
Turbine specific properties such as cut-in and rated speeds however, mean that 
the curve does not necessarily follow a quadratic and therefore the coefficient is 
not constant but itself varies as a function of flow speed. 

It is important to note, that the turbine properties derived in e.g. a CFD 
model, or from lab experiments, typically consider the placing of a single turbine 
in uniform background flow. Speed dependent properties are then expressed in 
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terms of the background velocity, which, because the velocity is slowed down in 
the presence of a turbine, is available as the undisturbed upstream velocity. In 
a finite width channel, blockage effects may also affect the resulting thrust curve 
but can be corrected for (see e.g. Garrett and Cummins [5], Whelan et al. [251 ) 
to derive the thrust curve for an idealised free-standing turbine. In addition, 
the results may be dependent on the turbulent properties of the flow. 

An approach followed in many models is to implement the thrust in the form 
of an equivalent drag force term. For depth-averaged models this effectively 
comes down to an increased bottom drag na [23i m 0 nu Three-dimensional 
models may implement the drag as a force over the entire water column 0 , or 
if the vertical resolution allows it the drag can be applied over a vertical cross 
section (e.g. Roc et al. [2D]), i.e. an idealised actuator disc. 

Since the thrust force is given as a function of the upstream velocity, it is 
important to consider what velocity to use for the equivalent drag force in the 
model. One option is to probe the numerical velocity solution somewhat up¬ 
stream of the turbine location. This however brings with it various difficulties 
such as the question of how far upstream is appropriate, or the fact that the 
flow upstream might not actually return to the uniform background flow con¬ 
dition that was considered in the CFD model, due to bathymetric changes or 
the presence of other turbines. Additionally, the use of a non-local velocity is 
not desirable for numerical and computational purposes: it makes it hard to 
treat the term implicitly (in the time-integration sense), potentially leading to 
time step restrictions for stability, and memory access outside of a fixed nu¬ 
merical stencil, or across sub-domains in domain-decomposed parallel model, is 
computationally inefficient. 

When enough mesh resolution is available, both in the horizontal and vertical 
dimensions, to resolve the flow through the turbine the relationship between the 
upstream velocity and the turbine velocity can be predicted using Linear Mo¬ 
mentum Actuator Disc Theory (LMADT). Using this relationship the quadratic 
drag law can be reformulated into a function of the local velocity, thus over¬ 
coming the difficulties and ambiguities mentioned above. This is the approach 
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followed in Roc et al. [20]. The typical width of a tidal turbine, order 20m, can 
however be orders of magnitude smaller than the spatial scales of the tidal flow 
so that resolving an individual turbine may become prohibitively expensive even 
in large-scale unstructured mesh models that allow for the efficient focusing of 
mesh resolution. Also this approach requires the alignment of the mesh with 
the position and direction of the turbine, thus limiting the flexibility to quickly 
evaluate different turbine positionings and angles. 

If the mesh resolution available is such that computational cells are much 
larger than the turbine scale, the drag force is necessarily applied over a larger 
area. In a typical implementation a constant drag is applied over a single cell 
(the cell that contains the turbine). If the cell size is in fact large enough it 
may be expected (this will be further investigated in this paper), that the local 
velocity is not actually affected greatly by the presence of the drag term since 
the drag force is “smeared” out over a large area and the local cell velocity 
represents an average of the velocity in a large area around the turbine. In that 
case the difference between the undisturbed background flow and the local cell 
velocity may be neglected and the turbine can be implemented using a simple 
quadratic drag law, function of the local velocity. 

As will be shown in this paper however, when the mesh resolution is refined 
such that mesh distances are closer to the turbine scale, this approximation is 
no longer tenable as the difference between upstream and local velocity becomes 
too large. As long as individual turbines are not resolved however, the approach 
in Roc et al. (20] is also not valid as the local velocity is still larger than the 
theoretical turbine velocity predicted by linear momentum actuator disc theory. 
In particular, for depth-averaged models the local velocity will remain higher 
than the actual turbine velocity even when the horizontal scales are sufficiently 
resolved. This is due to the fact that the drag acts on the entire water column 
and thus the depth-averaged model velocity will represent an average of the ac¬ 
tual turbine velocity and a higher by-pass velocity above and below the turbine. 
Even in three-dimensional models the drag force is often applied over the entire 
water column [5], or limited to one or only a few layers mm, and does not 
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necessarily give an accurate representation of the actual turbine cross-section 
and thus the model velocity where the drag is applied is not necessarily equal 
to the real turbine velocity. 

Here we demonstrate how the actuator disc computation may be modified 
to include the fact that the drag force numerically is applied over a different 
cross section than the actual turbine. Thus again an analytical relationship can 
be derived between the undisturbed upstream flow and the local cell velocity, 
and similarly the drag force can be reformulated as a drag law dependent on 
the local cell velocity. Like the approach in Roc et al. [20], this leads to a 
correction to the drag law, which in this case depends on the local cell width, 
but that nonetheless can easily be implemented in existing models, as will be 
demonstrated here for the Fluidity and MIKE 21 models. 

An alternative method for the parameterisation of turbines in large-scale hy¬ 
drodynamic models that also makes extensive use of actuator disc theory, is the 
line momentum sink method Hum- Actuator disc theory is used to express the 
effect of turbines, and more specifically an entire fence of turbines, as a relative 
head loss across the whole near-field flow pattern starting from the assumed 
uniform upstream flow at one end , down to the end of individual turbine wakes 
at the point where uniform flow is again achieved (within the far-field wake of 
the fence). This head loss is then applied as a jump condition across an edge, 
or multiple aligned edges within the computational grid using a Discontinuous 
Galerkin discretisation of the depth-averaged shallow water equations. The ad¬ 
vantage of this method that it incorporates a detailed LMADT treatment of 
the near-field effects, including blockage effects for multiple turbines in a fence. 
It does require however that these effects are treated at the sub-grid level, and 
is therefore only appropriate for hydrodynamic models with grid sizes larger 
than the length scale of the near-field/turbine wake (typically 10-20 turbine 
diameters) 17]. 

For any numerical modelling study it is important to look at the effect of 
changing the grid resolution on the results of interest. In the modelling guide¬ 
lines for tidal resource assessments in 14|, a range of grid resolutions is rec- 


5 


ommended depending on the stage of the resource assessment, ranging from 
kilometre scale for regional studies, down to a range of 500 m to 50 m for spe¬ 
cific site feasibility studies. Since the wake of a turbine is a three-dimensional 
phenomenon, it is not expected that an accurate description of the near-field 
flow can be obtained with a depth-averaged model. Nevertheless, such models 
should be capable of studying far-field effects. This relies on the correct forces 
and their effect on the large-scale flow being modelled correctly. As this paper 
shows however, the results of the standard enhanced bottom drag parameter- 
isation of the turbine thrust force will deteriorate as the mesh resolution falls 
below that of the near-field/wake length scale (« 200 — 300 m for a typical tur¬ 
bine) . The correction proposed in this paper ensures that consistent results can 
be obtained with grid resolutions smaller than the length scale of the turbine 
wake, all the way down to the turbine scale. 

2. Enhanced bottom drag formulation 

In this section we will describe the enhanced bottom drag parameterisation 
of turbines used in many models na 0 M 123 and demonstrate some issues 
with mesh dependency. We will do this within the framework of MIKE 21 [23], 
a depth-averaged hydrodynamics model widely used in the marine renewable 
industry, and an equivalent drag-based implementation in Fluidity, an open 
source, finite element modelling package dm12]. By comparing results between 
the two models we verify that the implementation in the closed source model 
MIKE 21 is indeed based on the same theory that underlies our implementation 
in Fluidity, and that the same issues are observed. 

The aim of the turbine parameterisation is to represent the drag force of the 
turbine on the flow, which is typically given as: 

F(u) = \pC t {\u\)A t \u\u, (1) 

here u is the flow velocity, p the density of sea water, C t the dimensionless 
drag or thrust coefficient, and A t the effective cross-sectional area of the tur¬ 
bine in the flow. The drag coefficient Ct may itself be a function of speed 
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due to turbine properties such as rating, pitch control and the use of a cut- 
in speed. As discussed in the introduction the drag law, often derived from a 
small-scale three-dimensional CFD model, is typically expressed as a function 
of the undisturbed background flow velocity, which in the case of an idealised 
domain corresponds to the uniform velocity upstream of the turbine. 

The depth-integrated shallow water equations (in conservation form) are 
given by 



( 2 ) 


(3) 


where H is the total water depth between bottom and free surface, elevated at 
a level z = r], u is the depth-averaged velocity, g the gravitational acceleration 
and Cb is the bottom friction coefficient. 

A local momentum balance in a fixed local horizontal area A is derived by 
integrating © over this area, multiplied by p\ 



(4) 


The second term represents momentum flux through the boundary dA. The 
third term can be rewritten as an integral of hydrostatic pressure around the 
three-dimensional water column below A. The last term represents a momentum 
sink term due to bottom friction. 

To implement the turbine thrust force through an enhanced bottom friction, 
c b —> Cb + Ct-, we need the additional momentum sink to be equal to the force, 
F(u) in ([I]). To address the question of which velocity u is used to compute 
F(u), in a first attempt we simply employ the local, depth-averaged velocity 
and average the force over the area A. Thus, we require that 



Combined with ([I]) , it readily follows that the enhanced bottom drag coefficient 
c t in this case should be set to: 



( 6 ) 
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Since we consider the parameterisation of turbines in hydrodynamic models 
where mesh distances are larger than the size of an individual turbine, the force 
is applied over the smallest area possible, typically the area of a single mesh cell. 
Thus the area A in ([6]) corresponds to the cell area over which the enhanced 
drag coefficient is applied. In models where the cell area is much larger than 
the turbine cross section A t , the additional drag is small and therefore the 
presence of the turbine will not have a large effect on the numerical solution for 
u in that cell. As an example, for typical values of Ct = 0.6, a mesh distance 
Ax = 200 m and turbine diameter D = 18 m, if the drag is applied over a single 
square computational cell of Ax x Ax, we get 

Ct{U)= Ct 2^ ^°- 00122 ’ ( 7 ) 

which is only half of a typical value of Cf, = 0.0025 for the background bottom 
friction coefficient. 

Since the effect of the additional drag is relatively small it is to be expected 
that the assumption that the local velocity within the cell is close to the undis¬ 
turbed background flow is valid for relatively coarse resolution models, and can 
therefore be used in the averaged force in the right-hand side of <§• As the 
resolution is increased however and the mesh distances become closer to the 
turbine scale, the drag is applied over a smaller area and the reduction in local 
flow speed may become much larger. Because of the quadratic dependency of 
the drag force on the flow speed, this may have a significant impact on the force 
that is applied in the model. 

3. Local velocity drop in idealised channel 

We investigate the mesh-dependent reduction in local flow speed in more de¬ 
tail in the following idealised set up: a turbine is placed in a rectangular channel 
of length 10 km and width 1 km. The depth at rest is set to 25m and a bottom 
friction of Cb = 0.0025, equivalent to a Chezy coefficient of 62.6 m 1 / 2 s _1 , is 
applied. At the upstream boundary a uniform velocity of 3.0 ms^ 1 is enforced. 



At the downstream end a Flather boundary condition is applied. The steady 
state solution without a turbine can be described as a balance between the free 
surface gradient and the bottom friction. The necessary free surface slope leads 
to a water level that is approximately 0.9 m higher at the upstream boundary 
than at the downstream boundary. The decrease in water depth H along the 
channel, in combination with the continuity equation, leads to an acceleration 
along the channel, with the speed increasing from 3.0 ms -1 to ss 3.12 ms” 1 
downstream. In a separate computation of the hydrodynamics without a tur¬ 
bine, it was verified in both Fluidity and MIKE 21, that the background flow 
velocity at the turbine location, halfway the channel, is approximately 3.055 
ms” 1 . 

For the simulations with a turbine, the following turbine parameters were 
chosen: the thrust coefficient Ct = 0.6 with a turbine diameter of D = 16 m 
giving a turbine cross-sectional area of A t = 201 m 2 . The simulations were 
performed using both MIKE 21 and Fluidity on a series of identical triangular 
meshes with uniform resolutions starting at a mesh size of Ax = 320 m, dou¬ 
bling the resolution each time with the mesh size decreasing down to Ax = 20 
m. One extra, fine resolution mesh with a mesh size equal to the turbine diam¬ 
eter was then run, Ax = D = 16 m. For the parameterisation of the turbine in 
Fluidity the enhanced bottom drag approach described in the previous section 
was chosen. Although Fluidity here uses a finite element scheme with discon¬ 
tinuous piecewise linear velocity and piecewise quadratic pressure solution (the 
mixed P1dg~ P2 velocity-pressure element pair, see Cotter et al. 0) a piece- 
wise constant drag field was used to simplify the computations and to remain 
close to the numerics of MIKE which uses a finite volume scheme with higher 
order flux reconstructions. Although the exact details of the implementation in 
MIKE were not available, the results between Fluidity and MIKE were found 
to be close enough to extend the analysis based on the parameterisation used 
in Fluidity to that in MIKE 21. 

Figure [T] displays the obtained velocity in the cell in which the drag has 
been enhanced to parameterise the effect of a turbine. MIKE 21 employs a cell 
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Figure 1: The speed at the turbine location, inside the enhanced bottom drag cell, decreases 
with increasing resolution both in the Fluidity and MIKE results. 

centred scheme, so for this model we report the value in the centre of the cell. In 
Fluidity’s numerical scheme, P1 dg — P 2, the velocity is represented by a linear 
function in each cell which is discontinuous between the cells. Here, and in the 
rest of the paper, the presented results for the local cell velocity are obtained by 
taking the cell average. From figure[l]it can be seen that the obtained velocity is 
indeed highly mesh-dependent, and drops with increased mesh resolution. Since 
the square of this velocity is used to implement the drag term, a 10% drop in 
the local velocity leads to an approximate 20% drop in the drag force. 

In a model of a fully resolved turbine the local velocity is expected to drop. 
After all, the velocity through the turbine is known to be smaller due to mo¬ 
mentum exchange with the turbine, whereas the bypass flow around the turbine 
is expected to accelerate. The deceleration of the flow through the turbine can 
be estimated using linear momentum actuator disc theory (LMADT, see [S] for 
an application of this theory to tidal turbines). The theory assumes inviscid 
flow and a uniform upstream velocity uq. Furthermore, it defines a velocity 
ui through the turbine, and velocities 113 and 114 in respectively the wake and 
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bypass flow (see figure [2]). It also defines pressures: po for the upstream pres¬ 
sure, pi and P 2 directly on either side of the turbine, and a uniform pressure P 4 
downstream where the velocities U 3 and 1(4 are defined. At the same downstream 
location, the cross-sectional area of the wake flow is defined as A 3 . In addition, 
it defines the known cross sections A c for the total channel cross section and A t 
for the turbine cross section. 

Through selective application of the continuity equation, momentum con¬ 
servation and Bernoulli’s principle, seven equations can be derived for the un¬ 
knowns ui, U3, U 4 ,pi,P 2 ,Po and A 3 , given uq and P 4 as upstream and down¬ 


stream boundary conditions respectively (see Appendix A). These equations 
can be simplified greatly by assuming A t -C A c , which means no blockage ef¬ 
fects are taken into account. For this case, 114 = uq,P 4 = po and the velocity 


through the turbine can be computed as (cf. equation (A.22) in the appendix): 

u l = \(l + y/l-Ct'j U 0 . (8) 


For our idealised channel case considered above, we may compute u\ = 
2.49 ms -1 . As we will see however in the next section, the difference between 
upstream and turbine velocity is smaller in the numerical results because the 
drag force is spread out over an area with a larger width. As already discussed, 
this means that at very coarse mesh resolutions the velocity in the drag cell is 
hardly different from the upstream velocity. As the mesh resolution is refined 
however, and the force can thus be applied over a smaller width, this velocity 
will drop. This decrease in local velocity will continue with increasing mesh 
resolution, and only when the resolution is sufficient that the drag force can be 
applied numerically over exactly the same cross section as that of the turbine, 
e.g. in a three-dimensional model, should we expect this velocity to have reached 
the value of Ui computed from standard actuator disc theory. 


4. Predicting the reduced velocity in the enhanced drag cell 

To simplify matters, we first consider the case where the turbine is repre¬ 
sented by a square area of enhanced drag instead of a triangle. Additionally, 
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Figure 2: Approximation of the enhanced drag formulation by actuator disc theory. An 
upstream velocity uq is assumed to reduce to a “turbine” velocity u ] inside the square in which 
the enhanced drag is applied. The effect of the enhanced drag is assumed to be equivalent 
to an actuator disc of width A y, the width of the cell in the direction transverse to the 
flow. The relationship between uj and uq can be estimated using actuator disc theory which 
involves eliminating wake and bypass velocities U 3 and 114 from a set of algebraic equations 
derived from selectively applying mass and momentum conservation and Bernoulli principles 
(see |Appendix A[ ). 

we assume that this square area is aligned with the flow. To this end we create 
a series of meshes with the same resolutions Aa; = 320 m to Ax = 16 m as in 
the previous section, but with an embedded square centred around the turbine 
location, of dimensions Ax x Aa:. The square is divided into two triangles, and 
outside the square an unstructured triangular mesh of approximately uniform 
resolution is created, as indicated in figure [2j 

When we neglect variations in the streamwise-direction (here denoted as 
the axdirection), the model results should correspond to those of an infinitely 
thin actuator disc as is considered in actuator disc theory. The actuator disc 
modelled by this shallow water model has a cross-sectional area of AyH. Here, 
and in the rest of the paper, Ay is the width of the drag area, in the cross¬ 
stream direction. In this section in particular Ay = Ax. Since we consider 
mesh resolutions where Ay > D 1 and additionally H > D, this “numerical” 
cross section will be much larger than the actual cross section A t . Therefore, 
to predict the results of the shallow water model using actuator disc theory, 
we should be careful to use the cross section AyH applicable to this model. If 
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we neglect any variations within the horizontal square, we may then hope to 
predict the velocity within the square as the disc velocity u\ from this modified 
actuator disc theory calculation. 

Following the assumption made above ([5]), the magnitude of the force applied 
in the enhanced bottom drag approximation is given by: 

F = \pA t C t u\. (9) 


Note that here we need to use the actual turbine cross section A t as that is the 
user input in this formulation to calculate the enhanced drag Ct in ([6]). Further 
we assume that the velocity that is used to compute the force in this approx¬ 
imation, which is simply the local velocity in the drag cell, will be accurately 
predicted as the velocity u\ in the modified actuator disc theory that follows 
below. 

Following the steps in the derivation of (Js|) , (A.221 in the appendix, but now 
applied to an actuator disc of cross section A t = AyH, we first define a modified 


c t . 


( 10 ) 


thrust coefficient (cf. (A.16) in the appendix): 

£ , = F = A t 

\pA t v? 0 A t u 2 0 

Following the same derivation of ([8]) , we then obtain a relationship between the 
local model velocity iti and the upstream velocity uq if in Q we replace Ct 
with Ct- This gives an expression for the ratio u\/uq than can be substituted 


in (101, to give: 


C t = f 


1 +VI -Ct 


Ct. 


After some algebraic manipulation^] this can be reworked to 


C t = 


At 

T t Ct 




in) 


( 12 ) 


1 


the authors made use of SymPy, a python library for symbolic mathematics: www.sympy. 


org 
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Finally, the relationship between the local velocity u\ within the cell that the 
enhanced drag is applied in, and the upstream velocity uq is given by 


u i = 


1 4- i At C 
1 + 4 A t C > 


-U 0 . 


(13) 


Figure [3] shows that the speed predicted by (131 closely follows that com¬ 
puted with Fluidity. Note that the results here differ from the Fluidity results 
in figure [l] This is because in figure |T| the drag is applied over an arbitrary tri¬ 
angle in an unstructured, triangular mesh generated by the mesh generator with 
a characteristic edge length set to the value of Ax on the x-axis. The meshes 
used for the results here, figure [l] are also unstructured, triangular and use the 
same characteristic edge lengths, but incorporate a square, consisting of two 
triangles, with dimensions Ax x Ax over which the drag is applied. Comparing 
the two figures, it can be observed that the drag being applied over a square 
area, aligned with the flow direction, leads to a different relationship between 
the upstream velocity and the velocity within the drag area, than when the 
drag is applied over a triangle. For this reason, in the following we will derive 
two different corrections to the enhanced bottom drag formulation for these two 
cases. 


5. Turbine correction for square cells 

We have shown that actuator disc theory, using the width of a square en¬ 
hanced drag cell and water depth, can accurately predict the relationship be¬ 
tween upstream and local cell velocities. This can be used to reformulate the 
drag force applied in the cell to be a function not of the local cell velocity, but 
effectively of the upstream velocity. Instead of applying the force in ([ 9 ]), which 
is based on neglecting the difference between upstream and local velocity, we 
want to apply the force 

F = \pA t C t ul, (14) 

where uq is the upstream velocity which is not readily (and locally) available. 
This expression is the same as in standard actuator disc theory, except now we 
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Resolution Ax in m 


Figure 3: The speed inside a square drag cell decreasing with increasing resolution. Results 
are model outputs from Fluidity. The plotted speed is the average value over the square area. 
The decreasing cell speed can be accurately predicted using EH derived from actuator disc 
theory. 


need to take into account that this force is not applied over the cross-section 
At but over a cross-section A t = A yH. Thus we obtain a modified thrust 
coefficient: 


dt := = f Ct • 

2pAtU 0 At 


(15) 


Note that this modified thrust coefficient differs from the one in the previous 
section, used to predict the results in the unmodified enhanced drag formulation, 
as we now assume that the correct force is applied. 

Assuming U\ is an adequate estimate for the local velocity in the cell with 
enhanced drag c t and cell area A, the force applied by the enhanced drag is 
given by (cf. the left-hand side of ([5])): 

F = pActuf. (16) 

After updating ([8]) to use the modified thrust coefficient C t . we can substitute 
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it here to make F a function of the upstream velocity u$: 

2 

“o- 


F — pAc t t 1 + V 1 — Ci 


( 17 ) 


To obtain the appropriate value of c* we simply equate this expression with the 


desired force in (14). This leads to: 

C t A t 


c t = 


2A 


( 1 + / 




2 ’ 


(18) 


In comparison with § from the standard enhanced bottom drag formulation, we 
have obtained an additional factor that corrects for the fact that we are using the 
local cell velocity instead of the upstream velocity. For coarse resolution runs, we 
have A t jA t —> 0, and thus we fall back, as expected, to the unmodified enhanced 
drag formulation, since the cell velocity is close to the upstream velocity. As we 
have seen for finer resolutions, still coarser than the turbine scale, the difference 
between cell and upstream velocities becomes significant. 

The correction derived above can also be applied to three-dimensional sim¬ 
ulations with a resolved turbine, where the drag force is applied in three- 
dimensions over a vertical cross-sectional area (actuator disc) with A t = A t 
and therefore Ct = Ct- The correction factor then simplifies to exactly that 
given in Roc et al. |20j . For the unresolved case however, both in two and three 
dimensions, the correction derived here not only corrects for the difference be¬ 
tween upstream and turbine velocity, but also for the difference between the 
actual turbine cross-section and the cross-section over which the drag is applied 
numerically. 

Returning to our idealised channel case, in figure [4] it is shown how the force 
in the standard enhanced bottom drag formulation applied to a square decreases 
with increasing mesh resolution. It is to be noted that the relative drop in drag 
force is larger than the relative drop in speed, due to the quadratic dependency 


of the force on the speed. Adjusting the drag formulation according to (18), the 
applied force is not only more accurate at coarse resolution, but also remains 
much closer to that computed from the upstream velocity directly as the mesh 
resolution approaches the turbine scale. 
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Figure 4: In the standard enhanced drag formulation for tidal turbines, equation j6|, the 
applied force is a quadratic function of the local velocity in the drag cell (here, a square area 
of Ax x Ax). As the mesh resolution increases, the local velocity drops, and therefore the 
force that is applied within the model decreases. Using the correction in ( | 18| ) however, the 
same force can be maintained more or less independent of resolution. 
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6. Turbine correction for triangular cells 


We now return to the case where the enhanced drag formulation is applied to 
a single triangular cell, not necessarily aligned in any way with the flow. Again 
we may approximate the applied drag by an actuator disc spanning the width 
of the triangle. In this case however, if we thus collapse the applied drag force 
to a single line, the amount of drag varies along the disc. 

We assume here that the streamlines run parallel through the triangle and 
use a local coordinate system where x is in the streamwise direction and 0 < y < 
Ay in the transversal direction, where Ay is the largest width of the triangle. We 
may subdivide the triangle into a number of streamtubes of infinitesimal width 
d y, which can be considered as rectangles Ax x dy, whose length Ax = Ax{y) 
is a function of y. When approximating this situation with actuator disc type 
theory, we make the following assumptions: 

1. The drag in each streamtube, which in the numerical model is applied over 
a length A x(y), is collapsed in the ^-direction and applied at a single point 
along the streamline, representing an infinitesimal actuator disc with cross 
section Hdy. 

2. The results in each of the streamtubes are independent of one another. 
This means that we take no blockage effects into account and assume 
laminar flow. 

For simplicity we first consider a triangle that is oriented in such a way that 
it is at its widest at y = Ay, in other words its top edge is aligned with the 
streamline at y = Ay (see figure [5]). Furthermore, we have A x(y = 0) = 0 in 
the bottom vertex, and A x(y) varies linearly for 0 < y < Ay. Its area can be 
computed as A = |Ax(Ay)Ay. The function A x(y) is therefore given by: 
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Ax ( y ) = ~^y- ( 19 ) 

The force applied in each streamtube is given by 

dF = Ax{y)&ypc t ui(yf, (20) 
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where u\{y) is the velocity through the streamtube. Similar to (10), we apply 
actuator disc theory where we assume that this force is applied over a cross 
section Hdy and obtain a modified thrust coefficient: 

dF 2A x(y)dy u\ 


C t ■■= 


\pHdyul Hdy 


jc t 


( 21 ) 


Following the same steps as in equations (10)—(13) we may derive the following 


relation between u\{y) and the upstream velocity uq: 

1 1 


ui{y) = 


1 i 1 Ax(y)dy „ 
1 + 2 Hdy C t 


U 0 = 


1 


HAy 2 y 


( 22 ) 


The varying width A x(y) thus leads to a variation of the velocity ui(y) for 
0 < y < Ay. In the computer models the accuracy of this variation is limited 
by the numerical approximations employed. 

In MIKE, the underlying discretisation is based on a piecewise-constant ve¬ 
locity in each cell. To estimate the cell average obtained in the model we 


therefore evaluate (22) in the centroid at y = |A y, which gives: 


i,2 Ac t 
L w 3 HAy 


U 0 . 


(23) 


For the case where the triangle does not have one of its edges aligned with a 
streamline, we may consider splitting the triangle into two triangles that share 
an edge that is aligned along the streamline (see figure [5]). The length of this 
shared edge is the maximum width A;r max of the triangular drag cell in the 
streamwise direction. The area of either of the two triangles that the cell is 
split into, can be computed as ^ 1,2 = ^A:r max A yi t 2 , where Ay\ 2 is the height 
of either triangle. Therefore for each of the triangles we have A\^/Ayi^ = 


1 Ax 


Thus if we apply the same enhanced friction coefficient c t in both 


triangles, it follows that the estimate (23) for the cell average of u“ IKB is the 
same in both triangles: 

(24) 


1 


1 max 

1 “t" 3 H 


Uq. 


Moreover, if we define the overall cross-stream width of the original combined 
triangle as Ay = Ayi + Aj/ 2 , we again have A/Ay = ^Ax max . Thus, in the 
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Figure 5: Left figure: a triangle with its top edge aligned with the streamlines. A coordinate 
reference frame is chosen, with 0 < y < Ay the coordinate in the cross-stream direction. The 
width Ax of the triangle in the streamwise direction, varies as a function of y , starting at 
Ax(y) = 0 at y = 0, and reaching it maximum width Ax(y) = Ax mELX at y = Ay. Right 
figure: a non-aligned triangle can be divided in two triangles that share an edge that is aligned 
with the streamlines. In this case, the maximum width Ar ma x is the length of the shared 
edge. 


actual model where the original, non-aligned triangular drag cell is not split, we 


can use the same equation (23) for the estimated average velocity of the entire 
cell as we did for the aligned case. 

Using this estimated average, the force applied in the model is then: 


F = Apc t (ur KE Y = Apct 


1 


1 + 


2 Ac t 


Ur,. 


(25) 


3 HAy , 


By equating this to the desired force ([T]), we may derive a quadratic expression 
for c t 


- 2 A 2 A t C t c 2 t + A (9H 2 Ay 2 - 6 A t C t HAy ) c t -\A t C t H 2 Ay 2 = 0 (26) 

In Fluidity, the PlDG _ discretisation prescribes a linear variation for velocity. 


Thus we approximate (221 by evaluating it at y = 0 and y = Ay and assuming 
a linear variation in between: 


r (») = i - 


l + ^A„ 


U 0 


(27) 
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The force applied in the model can be found by integrating: 


/■Ay 

F= / Ax{y) pet (u^ luidity (y)) 2 dy (28) 

Jy =o 

= - 4 «( 1 -!( IT W) + i(rykr) )< (29 > 

Equating with the desired force in ([I]) this time results in a cubic expression for 
c t : 


A 3 c d t + {AHAy - 3 A t C t ) cf 

+ 6 A (H 2 Ay 2 - A t C t HAy ) c t - 3A t C t H 2 Ay 2 = 0. (30) 

In case the triangular drag cell does not have an edge that is aligned with the 
streamlines, we may again consider splitting it into two triangles with a shared 


edge that is aligned with the flow. Here however, (27) does not predict the same 
linear function for w^ luldlty in both triangles, since although A/Ay = ^Ax max 
is the same, the value for Ay in the denominator of y/Ay is different for both 
triangles, and due to the different orientation of the top triangle, the sign of the 
gradient of iq luldlty with respect to y will be opposite. The combined piecewise 
solution is therefore not supported by the underlying discretisation. However, 


we did find that when using the value of c t found by solving (30), the discrete 


model gave results that varied only slightly for different orientations of the 
triangular cell. 

The results in hgure[6]indicate that again the force applied in the unmodified 
enhanced drag implementation, in Fluidity and MIKE reduces significantly with 
increasing mesh resolution. A modification to the enhanced bottom drag c* was 


derived in this section, solving for c t in (26) and (30) for MIKE and Fluidity 


respectively, that is shown here to lead to a force that remains close to the 
desired value. The correction in MIKE was implemented by first finding the 
value for c* from (26) and then working back from ([6]) to compute what value 
of C t should be entered in the GUI to achieve this value in MIKE. 
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Figure 6: Results for the enhanced drag formulation with the drag applied in a single triangular 
cell as implemented in both Fluidity and MIKE 21. As in figure |4| which show the results 
for the square case, the force applied decreases significantly with increasing mesh resolution. 
Applying the correction for ct however, given by solving ( |26| (MIKE 21) or ( |3Q[ ) (Fluidity), 
the force can be kept more or less constant and much closer to the desired value. 
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7. Power production 


The correction to the enhanced drag formulation, derived in this paper, is to 
ensure that the correct amount of momentum is extracted from a shallow water 
model. This means that the force F applied by the enhanced drag in the drag 
cell (or region) is an accurate approximation of the real thrust exerted by the 
turbine on the flow. The amount of energy taken out of the flow within the cell 
is given by: 

Pcell = Fill, model; (31) 

where rti,model is the velocity in the enhanced drag cell. As we have seen however, 
in the case where turbines are not fully resolved this velocity will be larger than 
the real velocity ui^urbine that goes through the turbine (as predicted by actu¬ 
ator disc theory). Therefore, the real power production Pturbine = Pity turbine 
will be smaller than the amount of power P ce ii taken out of the model in the 
drag cell. 

This discrepancy can be explained from the fact that part of the mixing 
losses are not modelled explicitly within the model, but occurs at the sub-grid 
scale. Following the analysis of Vogel et al. l 233 , the total amount of power 
taken out of the flow can be split as follows: 

Ptotal — Pturbine V Pmixing? (32) 

where P m i x ing takes account of the mixing losses due to a.o. shear between the 
wake and bypass flows. The total power can be computed as [21]: 

Ptotal = Fu 0 . (33) 


Therefore, as long as the model applies an accurate representation of the thrust 
force P, using the correction presented in this paper, and an accurate value 
for the upstream velocity Uq , the total power extracted from the flow in the 
model will be accurate as well. The fact that the power P ce u extracted within 


the drag cell, according to (31), is larger than Pturbine means that the mixing 
loss that occurs in the model (outside the drag cell) must be smaller than the 
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real P m ixing predicted by actuator disc theory. Therefore part of the mixing 
loss occurs within the drag cell itself. Thus P ce ii accounts for both the power 
Pturbine taken out by the turbine itself and additional losses that happen at the 
sub-grid level. 

Vogel et al. [21] considers the case where the drag of an entire farm is smeared 
out over an enhanced drag region, with the assumption that all mixing losses 
actually occur within this region. In that case it may be assumed that the total 
power extraction in the model is a good approximation of the total power extrac¬ 
tion predicted by actuator disc theory, so that the available usefully extracted 
power can be computed as a fraction of that using the same theory. 

For the case, considered in this paper, where individual turbines are modelled 
but are not necessarily fully resolved, part of the mixing losses are modelled 
explicitly. As argued above however, using the power extracted from the flow 
by the turbine parameterisation still leads to an overprediction of the usefully 
extractable energy. It is to be noted that in a shallow water model, even if 
an individual turbine is resolved in the horizontal mesh, with a minimum mesh 
distance smaller or equal than the turbine diameter D , the effective cross-section 
At = AyH will still be larger than the actual turbine cross-section A t . This 
is because the actual cross-section does not span the entire depth of the water. 
Thus, the velocity at the turbine in the model should be interpreted as a depth- 
averaged velocity that averages between the velocity through the turbine, and 
the bypass velocity above and below the turbine. This velocity is therefore 
expected to be higher than the real turbine velocity itself, and therefore the 
power extraction by the depth-averaged turbine-parameterisation will always be 
an overprediction of the actual power available to the turbine. The difference 
between these power values roughly corresponds to vertical mixing losses that 
are not explicitly modelled in the depth-averaged model. In the next section we 
will explain how the relationship between the upstream velocity and the local 
velocity in the model, derived in this paper, can also be used to predict the 
usefully extractable energy, excluding mixing losses, more accurately. 
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8. Implementation details 


In this section we summarise, how the analysis derived in this paper can be 
practically applied in existing models, in order to ensure that the correct force 
is applied on the flow and an accurate estimate of the available turbine power 
can be made. 


8.1. Turbine drag applied over a rectangular area 

For models where the turbine parameterisation consists of an enhanced bot¬ 
tom drag applied over a fixed, rectangular area A (e.g. [23]), we ma y use the 
analysis presented in section[5] Where existing models typically make no distinc¬ 
tion between upstream and local turbine velocity, they calculate the enhanced 
drag coefficient as Ct = CtA t /2A. Such implementations can be improved using 


the correction given by (18). The extra factor at the end of (18) can easily be 
included by the user in either Ct or A t , without the need for code modification, 
if these are the input parameters to the model. 

An additional complexity arises if C t itself is not a constant. This occurs for 
example if a cut-in speed and/or rating are applied to the turbine. In this case, 
Ct is typically given as a function (thrust curve) of the upstream velocity uq. 
In the model however only the local velocity u± is available. Using the formula 


U 1 — 2 ( 1 + V 1 — Ct ) Uo, Ct — . Ct: 


(34) 


however, it is straight-forward to transform a lookup table that gives the thrust 
coefficient for different values of uo, into a lookup table that is a function of u \, 
by computing ui for the given values of uq as a pre-processing step. 

For the computation of the power available to the turbine, we may use 


(A.23). Here, again we use (34) to derive the upstream velocity uq from the 


local cell velocity u±. Combining these two equations, we derive: 

2 (l + a/1 - C t ) 3 


Pt 


turbine 


<C t pA t u 1 . 


(35) 


^1 + \J \ — Ct 

Again, in the case that Ct is not a constant, a lookup table may be used to 
obtain the correct value of Pturbine foe each value of u\. 
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8.2. Turbine pcirameterisation in an arbitrary triangular mesh 

For models such as MIKE 21 and Fluidity that employ triangular meshes 
and which implement turbines through an increased drag applied within a single 
triangle, the theory presented in section [6] can be applied. In triangular mesh 
models where the drag force is based on a cell-averaged velocity, the value for 


the enhanced drag coefficient can be found by solving (26) for c*. Models that 
use a linear interpolation of velocities stored in the vertices, such as Fluidity 


should use the value of c t found by solving (30). The same approach could also 
be followed to implement a turbine in a single drag cell in Telemac 2D, where 
its Finite Element modus is expected to behave in a similar manner as Fluidity, 
using a linear representation of the velocity within a cell. 

In models, like MIKE, where the applied drag force and the associated coef¬ 
ficient Ct are not explicitly prescribed, the same effect can be achieved by modi¬ 
fying the value of C t . This is done by assuming the implementation is equivalent 
to the standard enhanced bottom drag formulation according to equation ([6]). 
Indeed the results in figure [l] where the standard drag implementation of Fluid¬ 
ity is compared with results in MIKE show that this is true to at least a good 
approximation. By providing MIKE with a modified value of Ct 

2Act 


C 


t, modified 


At 


(36) 


we can therefore create the effect of applying a value of Ct obtained from (26) 


without modifying the code. Note, that in equation (26) we use the original 
value of Ct for the real turbine. 

For non-constant Ct that is given as a thrust curve, MIKE (and similar 
models) use the local cell velocity u± instead of the upstream velocity to look 
up the value of Ct- This can be corrected by converting the upstream values uq 
in a uq —> Ct look-up table into cell velocities u\ using equation (23). 

To compute the power that can be usefully extracted by the turbine we again 


use (A.23) this time combined with (23), giving: 


^turbine — \ ^1 + \J 1 — Ct'j CfA t ^1 + 1 Ay} 


(37) 
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For finite element models, such as Fluidity, that consider a linear variation 


of the velocity within the cell we can use (271 which predicts the relationship 
between the upstream velocity and the velocity in the cell as a function of y. 
By first taking an average of the finite element solution u^ luldlty within the drag 
cell in the streamwise direction (x-direction), we can then use this equation to 
estimate the upstream velocity uq. This estimate may in practice still vary in 
the cross-streamwise direction (y-direction), so we take the cell average of its 


cube to obtain an estimate for Uq in (A.23). Combining all this gives: 


-Pturbine — \ + \/l ~ Ctj p 


C t A t 


r Ay 


L 


/y=0 


A x(y) 


Ai(j) . Fluidity 
“l 

A x(y) 


dx 


1 - 


l+HAy y 
Act A y 


dy. (38) 


8.3. Support structure 

The drag exerted on the flow by the support structure, e.g. pylons or tripods, 
can typically also be parameterised as a force that depends quadratically on the 
upstream velocity uq: 

-^support = (39) 

where A s and C s are the cross-sectional area and the drag coefficient of the 
support structure. To include this drag in the form of an enhanced bottom 
drag coefficient, we have to deal with the same issue of expressing this force 
in terms of a local velocity u\. Although the theory derived so far can be 
straightforwardly applied to any force in this quadratic form, we cannot simply 
derive the enhanced bottom drag coefficient that represents the support drag, 
denoted by c a , independently of c* and add them up. This is because the local 
velocity u\ is the depth-averaged velocity that will be slowed down by the two 
sources of drag simultaneously. 

For the drag parameterisation applied over a square area, the correct value 


for u\ is obtained from (341 by using 

A t C t + A S C S 


C t = 


HAy 


(40) 
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We can then derive a combined enhanced drag coefficient 


c t = 


A t C t + A S C S 


2 ’ 


( 41 ) 


2A (1+ 

that represents both turbine and support drag (cf. equation ©)■ The effec¬ 


tively extracted power is still given by (35) using the combined value of C t from 


(40), but only using the values for the turbine itself for Ct and A t . 

For models where the drag is applied over a triangular cell, the relation 
between uq and U\ is expressed in terms of the actual enhanced bottom drag 
coefficient c*. Thus if we include both turbine and support drag in c*, we 


can maintain equations (23) and (27) for models with cell-wise constant, and 
piecewise linear velocities respectively. The actual combined value of ct can then 


be found from the quadratic and cubic equations (26) and (30) respectively, by 
replacing A t C t with A t C t + A S C S . Finally, the extracted power is still given by 


(37), using only the turbine values for Ct and A t , but using the combined value 
of c t . 


9. Conclusions 

In order to accurately estimate the resource available to tidal turbines and to 
assess their impact on the hydrodynamics, it is important to accurately represent 
the drag force exerted by the turbines on the flow. In depth-averaged, and more 
generally under-resolved hydrodynamic models, one should keep in mind that 
the local model velocity at the turbine is different from both the upstream 
and the actual velocity passing through the turbine. The relationship between 
them is dependent on the mesh resolution, and in the case of depth-averaging, 
the ratio between the actual turbine cross section and the flow cross section 
spanning the entire depth. Therefore, although the use of the local velocity for 
the implementation of the drag force is computationally attractive, it is required 
to take these relationships into account to avoid spurious and mesh-dependent 
results. In addition, a better understanding of the relation between local and 
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upstream velocity is necessary for an accurate estimate of the power available 
to the turbine. 

Here we have presented the theory for a single, isolated turbine, and demon¬ 
strated that a correction based on linear momentum actuator disc theory taking 
into account the actual numerical cross section that the force is applied over in 
the model, can be used to obtain results that are consistent over a range of 
grid scales. It was shown that the standard enhanced bottom drag formulation 
results in a drag force that decreases with decreasing grid lengths, in particular 
when the grid size falls below the length scale of the turbine wake (roughly 
10-20 turbine diameters). With the correction the applied force can be kept 
constant to a large degree, thus ensuring that the effect of the turbine on the 
large scale flow is correctly modelled. 

The analysis for single, isolated turbines may be sufficient for sparsely pop¬ 
ulated turbine sites which see little interaction between turbines. It is generally 
recognised however, that in order to achieve the maximum available energy from 
certain sites, one needs to consider turbine configurations that benefit from local 
and global blockage effects HH, e.g. fence structures. The analysis in this paper 
could be extended to include blockage effects. Here again one should make a 
distinction between the influence of blockage on the relation between upstream 
and turbine velocities in reality, and the influence of blockage on the relation 
between upstream and local velocities in the model, in particular taking into 
account the difference in effective cross sections between reality and the model. 

With more closely packed turbines the representation of turbine wake struc¬ 
tures and wake recovery also becomes much more important. In addition, the 
turbulence characteristics may have a great impact on the performance of the 
turbines. As mentioned in the introduction, depth-averaged models will not 
be sufficient to accurately model these three-dimensional near-held effects. In 
further work we would like to explore however, how well these effects can still 
be approximated in depth-averaged models, possibly through parameterisation 
and tuning of horizontal turbulence models. Nonetheless, we recognise that in 
general it may no longer be possible to simply extrapolate from the results of a 
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single isolated turbine, and it may be required to study the effects of combin¬ 
ing multiple turbines in detailed three-dimensional CFD calculations and lab 
experiments. 
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Appendix A. Linear Momentum Actuator Disc Theory 

In this appendix we briefly review the main steps in the derivation of the 
actuator disc theory used in tidal turbine calculations. This is so we can refer 
to the relevant equations when the modifications, that take into account the 
numerical implementation details of the enhanced bottom drag formulation, 
are derived in the main text. These results can be found in e.g. Garrett and 
Cummins 0, or Whelan et al. [2BI . 

We consider a channel of cross-sectional area A c in which a turbine is located 
with cross section A t . We assume a uniform flow across the channel upstream 
of the turbine with velocity u o, the flow through the turbine is u±. Further 
downstream we define U 3 to be the velocity in the wake, and M 4 the bypass 
velocity. Furthermore we assume that at the point down-stream where M 3 and 
M 4 are defined we have a uniform water level 774 . The water level upstream 
is denoted by 770 , and the water levels just upstream and downstream of the 
turbine, associated with the pressure drop across the turbine are denoted by 771 
and 772 . 

First we formulate the conservation of mass for the flow through the turbine 
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and in the bypass flow 


A t u i = A 3 u 3 , 

A c uq = A 3 u 3 + (A c — A 3 )ui 


(A.l) 

(A.2) 


where A 3 is the cross-sectional area of the wake at the location where u 3 is 
defined. Here we neglect the influence of the water level on the cross sections, 
so that the cross-sectional area of the bypass flow is given by A c — A 3 . Inclusion 
of the dependency of cross section on the water level is only significant for high 
Froude numbers, with details given in |26j . 

The force F exerted by the turbine on the flow (and vice-versa), can be 
related to a conservation of momentum principle in the entire channel, or to the 
pressure drop across the turbine: 


F = A c pu 2 0 - A 3 pu\ - (A c - A 3 )pul + pgA c (g 0 - 774 ), 
F = pgA t (r]i -772), 


(A.3) 

(A.4) 


where g is the gravitational acceleration. Finally, applying Bernoulli’s principle 
along streamlines: 1 ) from upstream, where u 3 is considered uniform, to just 
before the turbine, where water level 771 is defined; 2 ) from just after the turbine, 
where water level 772 is defined, to downstream where a uniform water level 774 
is defined; and 3) in the bypass flow from upstream to downstream. This yields 
three more equations: 


h u o +9Vo = \u\ + gg 1 , 
\u\ +gp 2 = \ul +5774, 
\ul+gg 0 = \u\ + ggi. 


(A.5) 

(A. 6 ) 

(A.7) 


Assuming boundary conditions for uq and 774 , and an expression for A as a 
function of Uq, we have seven equations for seven unknowns: u 3 ,u 3 , W 4 , 770 , 771 , 772 , 
and A 3 . 


31 


General solutions 


The Bernoulli equations (A.5) to (A.7) can be rewritten as expressions for 
water level differences: 


gv i - gv 2 = gvo - gm + \ (ul - u \), 
gvo - gm = \ (ul - ul ), 

and thus 

gin - gv2 = \ (u 2 - ul). 

We can therefore rewrite the two expressions ( |A.3[ ) and ( |A.4 1 as: 

F = A 3 p (u 2 a - uf) - \A c p (ul - ul) 

F = \A t p(ul-ul) 


(A. 8 ) 

(A.9) 

(A.10) 

(A.ll) 
(A.12) 


Equations (A.2 1 , (A.ll I and (|A.12[) give three equations for the three unknowns 


u 3 ,U 4 and A 3 . Substitution of A 3 (m — u 3 ) =A c (u 4 — uq) from (A.2), in (A.Ill 
eliminates A 3 : 

F = A 3 p (U 4 - u 3 ) (U 4 + u 3 ) - \ A c p (ul - ul) 

= A c p (u 4 - u 0 ) (U 4 + u 3 ) - \A c p (ul - ul) (A.13) 

= A c p (U 4 - Mo) (m 3 + \u 4 - \u 0 ) . 

We can rearrange (A.12) and ( A.13[) in the following manner, respectively: 


A 2 C (u4 ~ w 0 ) 2 M 3 = A 2 C (u4 - M 0 ) 2 (u 2 - 
A 2 C (m 4 - mO) 2 M3 = - \ A c p (m 4 - m 0 ) 2 ) . 


We introduce the additional definitions, 

Ct := F 


and e := 


A t 


(A.14) 
(A.15) 

(A.16) 


\A t pul' ~ A c 

Note that we do not have to assume that F is actually quadratic in uo, so that 
C t is not necessarily a constant; it may still be dependent on mq. With these we 


can derive the following quartic polynomial in m 4 from (A. 14) and (A.15): 

, 2 


1 


C t e- 


U4 

M 0 


- 1 


M 4 

M 0 


- 1 


x 0 


-C t )= 0 . 


(A.17) 
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Finally, by (A.12): 


u 3 = \Ju\- C t ul , (A.18) 

and A 3 can be derived by again substituting dA~2l ) in (|A.11[) but this time to 


eliminate A c (it 4 — uO), so that 

F = A-ip (u 4 - u 3 ) (u 3 + \u 4 - \u 0 ) 


which in combination with (A. 12), gives: 

2^4 H - 2^*^ 


^3 = ^ ^ A t . 

zt 3 + 1>U 4 - ljUo 


(A.19) 


(A.20) 


2 4 2 

Zero blockage limit 

From the above, it follows that in the limit e —>■ 0: u 4 —> uq and thus r) 4 —> rjQ. 


In this limit, (A. 18) becomes 


u 3 


Vl -Ct 


t. Uq, 


and combining (A.20) and (A.l): 

\u 4 + ±u 3 


u i = 


u 3 + \u 4 - \u 0 
The energy yield then becomes: 


u 3 —> ^ (l + \/l — Ct'j Uq. 


(A.21) 


(A.22) 


P = 


Fui —t \ (l + 1 — Ct) CtAtpug. (A.23) 


The maximum yield as a function of C t is obtained by: 

1 - | C t + ^T^Gt 


d 

dCt 


(l + y/1 - C t ) C t 


VT^Cl 


= 0 


(fc* — i) 2 = 1 - Ct 


Ct = X. 


Thus the maximum power (assuming no blockage) is 
16 

f’max = 7 ^ ■ \A t pul « 0.59 • \A t pu\. (Betz limit) 


(A.24) 
(A.25) 

(A.26) 
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